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We analyze the range of validity of Thomas Fermi theory for describing charge density modulations 
induced by external potentials in neutral graphene. We compare exact results obtained from a tight- 
binding calculation with those of linear response theory and the Thomas Fermi approximation. 
For experimentally interesting ranges of size and density amplitudes (electron densities less than 
~ 10 11 cm~ 2 , and spatial length scales below ~ 20nm), linear response is significantly more accurate 
than Thomas Fermi theory. 
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I. INTRODUCTION 

The realization of single flakes of graphene - atomi- 
cally thin layers of carbon atoms packed in a honeycomb 
lattice - has made possible the experimental study of two- 
dimensional massless Dirac fermions [lj, [2|, l3| ■ Graphene 
is a gapless semiconductor in which the conduction and 
valence bands touch at two points - Dirac points - in the 
Brillouin zone Near either of these points the elec- 
tronic states are described by a massless Dirac equation, 
with eigenstates which are spinors due to the two-point 
basis needed to describe the honeycomb lattice The 
effective spinors of the wavefunctions are either parallel 
or antiparallel to the momentum, so that the states are 
chiral. 

For undoped graphene there is one electron per car- 
bon atom, and the system ideally should be everywhere 
charge neutral. In practice this is known not to be the 
case. Recent imaging experiments @ have demonstrated 
the existence of electron and hole puddles of densities 
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point. The existence of these charge puddles could be re- 
lated to the existence of mechanical ripples also observed 
in graphene sheets @, H, Q, which can cause modula- 
tion of the electronic charge [13, [H| > or to unintentional 
charged impurities in the substrate (l2. 115, 141 which 
can also generate electron- hole puddles 15 Jig 17j . The 
spatial correlation length of these puddles is of the order 

Of 1071771. 

Local density inhomogeneities can also be induced in 
graphene using miniature gates. In this way graphene 
p-n junctions have been experimentally realized [la . Il9l . 
|20j. Recent advances in the quality of graphene have 
made possible the fabrication of ballistic circuits with 
electrically controlled p-n i unctions [2ll |22| . 

The physical properties of graphene with such elec- 
tronic inhomogeneities depend strongly on the size and 
amplitude of charge modulation induced by external po- 
tentials. It is therefore important to understand how the 
ground state charge in graphene is distributed in their 
presence. Large inhomogeneous graphene systems have 



been studied theoretically using the Thomas Fermi (TF) 
approximation, which, as we discuss below, treats the ki- 
netic energy in a local density approximation [23| . Rossi 
and Das Sarma used a TF approximation with Hartree 
and exchange effects included to study the ground state 
of neutral graphene in the presence of charged impurities 
[l5l |. A more rigorous quantum mechanical treatment of 
the kinetic energy is possible, but its use limits consid- 
erably the system sizes which in practice can be studied 

As we will show below, because of the crossing of the 
chiral electron and hole bands at the Dirac point, the 
TF approximation does not correctly capture the charge 
response of neutral graphene to an external potential in 
many interesting situations. The purpose of this work is 
to analyze the range of validity of the TF theory near 
the Dirac point. We use a microscopic tight-binding cal- 
culation to compute the response of neutral graphene to 
electrostatic potentials, and compare these exact results 
both with linear response and with the TF approxima- 
tion. We will demonstrate that for experimentally inter- 
esting 6] ranges of sizes and amplitudes (electron densi- 
ties ~ 10 n cm -2 and spatial correlations ~ 20nm), sim- 
ple linear response results match exact results quite well, 
while results of the TF approach are much poorer. The 
failure of the TF approximation is related to the non- 
local character of the density response, and we shall see 
that a kinetic energy functional which correctly captures 
the linear response of neutral graphene to external elec- 
trostatic perturbations has a highly non-local nature. 



II. THOMAS FERMI FUNCTIONAL FOR THE 
KINETIC ENERGY. 

A. Formal Considerations 

Following Hohenberg and Kohn [24j , the total energy of 
the noninteracting system, E, may be written in terms of 
a kinetic energy functional T[n(r)] of the electron density 
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ra(r), 



E[n(r)} = J T[n(r)]dr + J 



V{r)n(v)dT 



(1) 



Here V(r) is the one-body external potential in which the 
particles move, and the density is defined with respect to 
the density of electrons in neutral graphene. The effect of 
electron-electron interactions in a Hartree approximation 
will be considered below in Section III. 

The TF theory assumes that the functional T[n(r)] 
is a local function of the density, and the form of the 
functional is chosen such that for a uniform potential V 
the minimization of Eq. [1] recovers the kinetic energy of 
a homogeneous system. For the case of Dirac fermions, 
the Thomas-Fermi kinetic energy functional is 



T[n] 
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where Up is the Fermi velocity of the carriers near the 
Dirac point. The minimization with respect to the den- 
sity must be carried out subject to the normalization 
constraint 



n(r)dr = uq 



(3) 



where no is the average electron density measured rela- 
tive to that of undoped graphene, and S is the sample 
area. Minimizing Eq[T] yields the relation 
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where u = hvFy/nno is the Fermi energy of the cor- 
responding homogeneous system. Defining k max {v) — 
y / n\n TF (r)\, we find that the carriers have higher ki- 
netic energy where the potential energy is lower, and 
vice- versa. Eq. [3] van be viewed as the relation between 
the local maximum momentum k max (r) and the exter- 
nal potential V(r) obtained from a classical equation of 
motion. 

An interesting and important consistency check of the 
TF approximation [24[ is that the response function of 
the system should be directly related to the second func- 
tional derivative of the energy. For a non-interacting sys- 
tem this takes the form 
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where T indicates the Fourier transform, and XLin{q, no) 
is the wavevector dependent static Lindhard susceptibil- 
ity of the uniform system at density no. 

In graphene the Lindhard static susceptibility in the 
long wavelength limit has the form[14j, |25|, |26| 
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for q < y/nno, 

for n = 0. (6) 




FIG. 1: (Color onZme)Unit cell used in the calculations. The 
unit cell contains N atoms and the length of the unit cell is 
L — N/4ao. The external potential only depends on x and 
in this geometry atoms A and B with the same x coordinate 
have the same charge [27, . [28j] . ao is the lattice parameter of 
the triangular lattice. 



The TF kinetic energy functional, EqJ^l correctly re- 
covers the response function for doped graphene in the 
long wavelength limit, but it fails to describe the non- 
interacting compressibility at the neutrality point. In 
fact, the TF approximation predicts vanishing linear re- 
sponse at the Dirac point. This failure is in agree- 
ment with the general assumption of the TF theory that 
|Vn(r)|/[n(r)fc m(M: (r)] <c 1, which cannot be satisfied 
near charge neutrality. Moreover, as we discuss below, 
the response in the second of Eqs. [6] is inherently non- 
local, suggesting that the TF approximation must break 
down near charge neutrality. 



B. Numerical results. 

In order to quantify the effects of the failure of the TF 
approximation to correctly describe the linear response 
of undoped graphene, we numerically compute the elec- 
tron density of a net neutral graphene system in an ex- 
ternal potential, and compare the results with the TF 
approximation and with linear response results. We use 
a simple tight-binding Hamiltonian with nearest neigh- 
bor hopping, of the form 



(7) 



<hj> 



where Ci annihilates an electron at site Ri of the 
graphene lattice, t = ^j^p is the hopping paramenter, ao 
is the lattice parameter of the triangular lattice, and Vi 
represents the external potential at site R;. We perform 
the calculations in a unit cell illustrated in Fig. [TJ using 
periodic boundary conditions in both the x and y di- 
rections. The external potential and the induced charge 
depend only on the x coordinate. In the unit cell rep- 
resented in FigfT] atoms on both sublattices experience 
the same external potential, so there is no out-of-phase 
response from atoms on different sublattices [2^. 

We study the response of the system to the potential 



Vi = V cos{GXi 



(8) 
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FIG. 2: (Color online) Density profiles obtained with differ- 
ent approximations for perturbation amplitude Vo = 50meV 
and period 100ao. For exact calculations, t = 2.8eV and 
ao = 2.46A. Solid line is the exact result, dashed line indi- 
cates linear response result, and dash-dotted line is result of 
Thomas Fermi approximation. 



where Xi is the x-component of the position of the carbon 
atoms, and Vb is the amplitude of the perturbation. Fig[2] 
illustrates a typical result, the electron density induced 
by a potential of amplitude Vo = 50meV and period 
100a . Also plotted are the density as obtained in lin- 
ear response, n Lm (G) — XLin(G, 0)Vq, and from the TF 
approximation, Eq. 2] The density induced by this po- 
tential is of the same order as the densities of electron and 
hole puddles observed experimentally. Note that the lin- 
ear response reproduces the exact result rather faithfully, 
whereas for this potential the TF approximation under- 
estimates the response. Moreover, the TF approximation 
displays plateau-like features when passing through zero 
density, which are an artifact of the approximation [23|; 
they appear because TF theory grossly underestimates 
the ability of the system to screen when the local chem- 
ical potential is near the Dirac point. The plateaus may 
be understood more formally by substituting the pertur- 
bation Eq. [5] into Eq. 2] and expanding in harmonics, to 
obtain 
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The large cos 3Gx harmonic leads to the plateau-like be- 
havior when crossing the Dirac point. 

In Fig. [3] we compare the maximum electron density 
at x — 0, obtained both from the exact calculation, and 
in the two different approximation schemes, as a func- 
tion of Vo, for different periods of the external potential. 
For small periods and small Vo, the linear response re- 
sults follow the exact results rather closely. TF theory by 
contrast underestimates the response of the system. For 
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FIG. 3: (Color online) Maximum induced electron density 
(density at x—0) as function of the amplitude of the exter- 
nal potential. The external potential has the form V(x) = 
Vo cos Gx. (a), (b) and (c) correspond to values Gao = 7r/200, 
Gao — 7r/100 and Gao = n/50 respectively. Continuous lines 
are the exact results, dashed lines are the linear response re- 
sults, and dash-dotted lines are the results obtained in the 
Thomas Fermi approximation. In the calculations we use the 
values t = 2.8eV and a = 2.46A. 



small enough Vb and large G, linear response is able to 
properly capture the non-local nature of screening in this 
system. For large wavelengths and external potentials 
non-linear contributions to the response become impor- 
tant, and may be captured by the TF approximation in 
any average way (Fig. [3]Ja).) From the numerical results 
we estimate that, in the absence of electron-electron in- 
teractions, linear response is more reliable than TF when 
n Lm > n TF . For large perturbations the exact density 
response oscillates around the TF result. These oscilla- 
tions are induced by zero modes created by the external 
potential in graphene(30j. which cannot be captured by 
a local approach such as the TF approximation. 

For the charge density modulation amplitudes ob- 
served experimentally, ~ 10 n cm -2 , the length scale 
for which linear response is more reliable than the TF 
approximation is larger than the size of the observed 
electron- hole puddles [6|]. Furthermore, from the geom- 
etry of the multiple gated graphene devices in Refs. [U 
and I22I we find that the width of the depletion regions 
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in the p-n junctions [23j are also smaller than the length 
scale where linear response is applicable. More generally, 
our results indicate that for density modulations up to 
10 12 cm~ 2 on length scales up to 20nm, linear response 
results are significantly more accurate than those of the 
TF approximation. This conclusion agrees with results 
presented in Fig. 2 of Ref. [l|| where the authors find re- 
sults which are consistent at semi-quantitative level with 
a linear screening theory. 



III. HARTREE INTERACTION 

A. Formulation in Terms of Linear Response 

Any modulation of electric charge produces a change 
in the energy associated with the repulsion between elec- 
trons. If one is interested in the long- wavelength static re- 
sponse of the charge density to a potential inducing such 
a modulation, the most important effects of the electron- 
electron interaction can be captured by the Hartree en- 
ergy. This may be written in the form 



E H = 



2~e 
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where e is the average background dielectric constant. 
The strength of the Coulomb interaction is given by the 
dimensionless parameter 



(11) 
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For graphene on a conventional SiC>2 substrate, e w 2.5 
and a » 0.9. For substrates with larger e such as HfC>2 
or liquid water, the values of a can be much smaller. We 
note that in principle one may improve upon the Hartree 
approximation by including exchange correlation effects, 
but for chiral Dirac fermions these appear to be rather 
small 0. 

Since we will consider perturbations with amplitudes 
and periods such that the exact non-interacting result 
coincides nearly perfectly with that of linear response, 
we expect that the inclusion of the Hartree term leads 
only to linear screening of the external potential. In this 
case the induced charge coincides with that obtained in 
the Random Phase Approximation (RPA). In reciprocal 
space this means 
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is the two dimensional Fourier transform 



of the Coulomb interaction. 

In the TF approximation the electron density is ob- 
tained by minimizing the kinetic functional, Eq. [T] to- 
gether with the Hartree energy Eq. \W\ with respect to 
the density. In Fig. 0] we compare the spatial maxi- 
mum electron density obtained in the RPA to the TF 
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FIG. 4: (Color online) Maximum induced electron density 
as a function of the amplitude of an external potential of pe- 
riod lOOao, in the Hartree approximation, (a) corresponds 
to a weak electron-electron interaction, a— 0.5 and (b) to 
a stronger one, a—1. Dashed lines represent RPA results; 
dash-dotted lines are results obtained by minimizing an en- 
ergy functional containing both the TF approximation to the 
kinetic energy and the Hartree form of the Coulomb interac- 
tion. 



approximation as a function of the amplitude of the ex- 
ternal potential. The Hartree interaction screens the ex- 
ternal potential, so that the induced charge density de- 
creases with increasing electron-electron interaction pa- 
rameter a. As in the non-interacting case we see that 
the TF approximation underestimates the response at 
small Vb- For physically relevant values of a, we see that 
the TF approximation is not quantitatively reliable in 
describing the response of neutral graphene to external 
potentials that generate density fluctuations of magni- 
tude 10 12 cto~ 2 or below, within length scales of about 
20nm. 



B. Electric fields in a p-n junction 

Ballistic transport in graphene p-n junctions is due to 
Klein tunneling of the massless electrons. Cheianov and 
Falko [3l| showed that the ballistic resistance per unit 

width of a graphene p-n junction is R — | where 
E is the assumed uniform electric field at the junction. 
Note the resistance decreases as the electric field at the 
interface decreases. This electric field depends on the 
screening properties of graphene near the Dirac point. 
Zhang and Fogler [23] proposed that the electric field 
in the depletion region separating the electron and hole 
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FIG. 5: (Color online) Electric field as a function of position 
for a periodic external potential of amplitude Vo=0.1eV and 
period 100ao. Coupling constant is taken to be a = 0.5. 
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regions is enhanced due to the limited screening capacity 
of Dirac quasiparticles. 

In order to study the difference in computed values 
of F, the electric field in the depletion region, using the 
TF approximation and linear response theory, we have 
calculated the electric field for a cosine-shaped external 
potential EqJS] This potential creates periodic electron 
and hole regions separated by p-n interfaces. In Fig[S] 
we plot the electric field as a function of position, as ob- 
tained in the TF approximation, and in the RPA (the 
latter being essentially an exact solution of the Hartree 
approximation.) For comparison we also plot the applied 
electric field, E ext = —GVq sinGx. The results presented 
are for a=0.5. The p-n and n-p interfaces are located at 
x—25ao and x=75ao, respectively. At these points the 
values of the electric field are maximal. In the linear 
calculation the electric field F can be calculated analyt- 
ically, yielding the result F = Vq G/(l + 7r/2a), so that 
the external electric field is reduced by a factor (l+n/2a) 
by the screening. In the TF approximation a numerical 
minimization is required to obtain F. In the range of 
validity of the linear approximation we find that the TF 
approach predicts much weaker screening of the external 
field than the RPA. 

In Figs. Oa) and^b) we plot the values of the electric 
field at the p-n junctions, normalized to the external field, 
as a function of the applied electric field, for two different 
values of a. The screened electric field at the interface 
obtained from the TF theory is larger than that obtained 
in the linear response theory (RPA), as expected from 
the above results. We see that the TF approximation 
significantly overestimates the total electric field at the 
p-n junction. 



FIG. 6: (Color online) Electric field at the p-n junction 
centers as a function of the external field, (a) corresponds 
to a weak electron-electron interaction, q=0.5 and (b) to a 
stronger one, a=l. Continuous lines corresponds to the RPA 
results. Dash-dotted lines are the results obtained in the self- 
consistent Thomas-Fermi approximation. 

IV. SUMMARY AND OBSERVATIONS 

The Thomas Fermi approximation is relatively inaccu- 
rate for describing density modulations for wavevectors 
that are not too small, and external potentials which are 
not too large, in undoped graphene. Quantitatively this 
region of failure of the TF approximation appears to ap- 
ply to the observed density fluctuations of the electron- 
hole puddles that appear in the single electron transistors 
spectroscopy. It also appears to be problematic for esti- 
mating the electric field in a graphene p — n junction. 
The reason for its failure is its inability to capture the 
intrinsically non-local response of neutral graphene. We 
find that the application of linear response theory (RPA) 
in this regime is far more quantitative. 

It is interesting to note that one may adopt a non-local 
kinetic energy functional to produce a correct result for 
Eq. [5l This takes the form [32[ 

T««->(r)] = ^L[ d r(dr> (13) 
tv J J |r — r'| 

This kinetic energy functional is formally the same as the 
Hartree form of the interaction energy, highlighting the 
marginal nature of 1 jr Coulomb interactions in undoped 
graphene [4[ . Its long-range nature strongly suggests the 
difficulties of a local approximation such as TF that we 
find. 

To improve upon the TF approximation one can for- 
mally compute first order gradient corrections to the den- 
sity using a WKB approximation applied to the Green's 



G 



function [33|.The result ,34], however, has singular be- 
havior near zero momentum, and moreover depends lo- 
cally on both the density and its gradient, and so cannot 
produce corrections where the density is maximum and 
where TF has significant errors. 

Finally we note that Eq. [13] may be used to develop 
a criterion for which one expects the TF approximation 
to fail. Using the result of the linear response density in 
Eq. [13] gives an estimate for the energy density expected 
from the non-local contribution to the energy. Compar- 
ing this to the TF energy density (Eqs. [5] and 01, we 
expect to the latter to be larger if the TF approximation 
is to be valid. This yields the criterion Vb/G > r/hvF, 
where 77 is a geometric factor of order 1, for which the 
TF kinetic energy dominates over non-local contributions 
to the energy. Notice this means that, for fixed length 



scale 1/G, the TF approximation will always fail for suf- 
ficiently small potential scales Vq. 
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